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Abstract 

Using a simple and generic molecular dynamics model, we study the 
damage in a disc of interacting particles as the disc fragments upon im- 
pact with a wall. The damage, defined as the ratio of the number of bonds 
broken by the impact to the initial number of bonds, is found to increase 
logarithmically with the energy deposited in the system. This result im- 
plies a linear growth with damage for the total number of fragments and 
for the power law exponent of the fragment size distribution. 

PACS numbers: 46.50.+a, 62.20.mm 
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1 Introduction 

Because of its importance in natural phenomena and industrial processes, frag- 
mentation is of great interest. To a certain extent, fragmentation phenomena 
depend on the particular features of the object that is broken. However, some 
generic behaviours seem to be shared by fragmenting systems whatever their 
size, material, or typical interaction energy. For instance, as observed in a 
large variety of experiments [TJ [2j El [4j O [71 [HI El HH] and natural phenom- 
ena [TJ HU [T2] , the fragment size distribution frequently exhibits a power law 
behaviour, the origin of which is still unknown. 

Despite its simple and rather everyday aspect, fragmentation is a complex 
and distinctly far-from-equilibrium process which makes it difficult to under- 
stand theoretically. Attempts have been made to derive a power law distribu- 
tion using analytical models, such as sequential fragmentation [5] and energy- 
balance theory [T3]. However, these approaches do not succeed in explaining 
the variety of values of the power law exponent observed in experiments. On 
the other hand, fragmentation is beneficially studied with computer simulations 
[HI HH HH HZl HH US HOI HH HSl Indeed, this approach can grasp 

the complexity and the dynamical properties of fragmentation by taking into 
account the various parameters that may influence the process — and consider 
even them one by one. 

The main goal of our studies is to highlight and understand the behaviours 
observed in a wide range of different fragmenting systems. To encompass this 
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variety, we recently proposed an elementary molecular dynamics (MD) model 
[24| which, due to its simplicity, provides a very generic frame of reference for 
fragmentation studies. Therefore, this model is not representative of any specific 
material. We applied it to investigate the fragmentation of a two-dimensional 
disc of interacting particles upon impact with a wall [24 . One of the main 
results of this study is, that the power law exponent of the fragment size distri- 
bution increases logarithmically with the energy deposited in the system. This 
behaviour is in agreement with experimental results [3J [TU] and we expect it to 
be generic for fragmentation phenomena. 

In the present paper, we extend this study by addressing the issue of damage 
occurring in a fragmenting disc. In particular, we will show how the breaking of 
bonds at the microscopic level is related to macroscopic quantities such as the 
fragment size distribution. Dynamical quantities dealing with the spreading of 
damage during the fragmentation process will also be discussed. 

The paper is organized as follows: in section 2, we describe the model and 
the simulation procedure ; in section 3 we review the main results concerning 
the fragment size distribution ; section 4 is devoted to the study of damage ; in 
section 5 we discuss the fragmentation energy. We conclude in section 6. 

2 Model and simulations 

In the present study, the fragmenting system is a disc made up of N — 1345 
particles placed on a two-dimensional triangular lattice. In the previous article 
[24j , we showed by a finite size scaling approach that the fragmentation features 
of the system are not sensitive to the number of particles. The cohesion of the 
system is ensured by a central two-body Lennard-Jones type of potential 



Here ry is the distance between particles i and j and the two constants, e and a, 
are the depth of the potential well and the diameter of the particles, respectively. 
The three parameters, vo = 107.37, a = 80 and b = 78, were chosen to obtain a 
strong repulsion at contact and a very short range of attraction (around 0.1 a in 
addition to the particle diameter). These are the features that one expects to be 
dominant within a brittle solid at mesoscopic scales. Note that the behaviour 
of the system under fragmentation is not sensitive to a particular choice of the 
parameters vq, a, and 6, that is, to the range of attraction [24j . 

The fragmentation process was studied by means of 2D molecular dynamics 
simulations utilizing the Verlet algorithm [25} 126] at constant energy. The time 
step is chosen as St = 0.0001 to, ensuring the conservation of total energy up to 
0.001 %. Here to = \A<7 2 /m is the unit of time and m — 1 the particle mass. In 
the following, times, lengths, velocities, and energies are expressed in the units 
of to, a, a/to, and e respectively. 

The disc of diameter L — 40er is constructed, rotated by a random angle and 
launched into a solid wall with a given impact velocity V. The particles of the 
disc interact with the wall through the repulsive part of the potential v(i"ij) (for 
details, see [24). For each value of the impact velocity, we performed 1000 runs 
while sampling the initial angle of rotation uniformly. Simulations were stopped 
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at t = 100 (i.e. after 10 6 MD time steps), when the fragmentation process had 
already reached a steady state [21] ■ 

During the simulation, fragments are defined as self-bound clusters of parti- 
cles [27J . In other words, two particles are bonded if their relative kinetic energy 
is lower than the absolute value of their interaction energy. Once the distance 
between two particles i and j becomes larger than the attractive range (~ 1.2a) 
of the potential, their potential energy is set from e to O.OOle. This value reduces 
the interaction between these two particles to an almost repulsive one and thus 
the possible recombination is prevented. 

In this article, we measure the energy deposited in the system with a dimen- 
sionless control parameter 77, defined as the ratio of the initial kinetic energy to 
the potential energy of the disc: 

mV 2 

11 = 01 r 

^ ppotOl 

Here e po to = — 2.87e is the initial potential energy per particle in the disc. 
When the kinetic energy is large enough compared to the cohesion of the disc, 
fragmentation occurs. 

The same numerical experiment — a disc of interacting particles fragment- 
ing upon impact with a wall — has been previously investigated by means of 
more complex models. For instance, Thornton and his collaborators [18) used 
primary particles interacting through autoadhesive and frictional forces that 
act as functions of contact area. The particles, set initially in random locations, 
are pulled together by applying a centripetal force. Hence in their model, the 
fragmenting agglomerate has an irregular shape and contains inherent struc- 
tural defects. Likewise, in the model proposed by Kun, Herrmann, and their 
collaborators [201 [22 [22] , the primary particles are rigid randomly shaped con- 
vex polygons. These polygons are connected with beams that take into account 
elastic, shear, and torque interactions. 

The asset of the present study is the use of a minimal fragmentation model 
since it provides a generic frame of reference to which one can compare the 
effects of different parameters. Accordingly, the system considered in this article 
is perfectly ordered and homogeneous, but disorder and heterogeneities could 
be easily taken into account to investigate their effects. 

3 Fragment size distribution 

For a given amount of energy deposited in the system, we calculated the frag- 
ment size distribution — the number n(s) of fragments made up of s particles 
— averaged over 1000 runs. In addition to this, we calculated the normalized 
mean size of the largest and the second largest fragment — S max i and S rnaX 2 — 
as well as the total number of fragments mo- These last quantities are plotted 
in Fig. [T]as functions of rj at the end of the simulation (t = 100). As we found 
in the previous study [21], S max i decreases with r], whereas S maX 2 is peaked 
at rjt — 0.19 (Vt — 1.05). When the amount of energy deposited in the system 
is larger than this threshold value rj t , the fragment size distribution exhibits a 
power law behaviour with an exponent t: n(s) ~ s~ r . 

To estimate the exponent r, we employed statistical procedures proposed by 
Clauset et al. [2E] , which are based on the method of maximum likelihood. Note 
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Figure 1: The mean size of the first and second largest fragments, S max i/N and 
Smaxv/N, and the total number of fragments tuq/N as functions of n at time 
t = 100. The full line corresponds to a logarithmic fit (see eq. §Bj))- Dotted lines 
are a guide for the eyes. The threshold energy rjt is indicated by an arrow. 



that, as 77 increases to values much larger than rjt, the power law region of 77,(5) 
narrows due to the finite size of the fragmenting disc, and the determination of 
r becomes difficult. Nevertheless, the power law fit is quite good for rj t < n < 4, 
(see the fragment size distribution plotted in the inset of Fig. [2]). 

Furthermore, as illustrated in Fig. [21 the exponent r of the power law 
exhibits a logarithmic increase with the energy deposited in the system: 

t = a T \n— + P T , (2) 
Vt 

where a T = 0.26 ± 0.01 and (3 T — 1.34 ± 0.01. When t is given as a function of 
the impact velocity, the coefficients a T and (3 T are slightly different than those 
obtained in the previous article [21]. This is due to the fact that, in the present 
paper, we use a more rigorous method [25] for fitting. 

Similar to the size of the largest fragment, the total number of fragments 
mo = Yl s n(s) is a measure of the degree of break-up in the fragmentation 
process. Figure Q] displays also mo as a function of 77. As expected, the total 
number of fragments increases with the initial kinetic energy. Moreover, mo 
presents a linear behaviour with r\ for 77 < 1, and a logarithmic behaviour for 
1 < 77 < 10, that is, up to the highest impact energies considered in this work: 

= a mo \nn + [3 mo . (3) 

Here a mo = 0.26±0.01, i.e., the same coefficient as for r, and [3 mo — 0.23±0.01. 
Note that Behera et al. [21] found the same linear and logarithmic behaviours 
while studying the impact fragmentation of a disc with a more sophisticated 
model, that takes into account also the elastic, shear and torque interactions 
between particles. 
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Figure 2: Power law exponent r as a function off). The full line is a logarithmic 
fit (see eq. (ify). Cumulative fragment size distribution for V = 3 (rj = 1.57) 
at t = 100 is plotted in the inset (the dotted line is a power law function with 
t = 1.87). 



4 Damage 

The fragment size distribution discussed in the previous section reflects the 
degree of damage caused by the impact. Furthermore, the increase of the to- 
tal number of fragments with the impact velocity would be, and is, a direct 
and intuitive representation of increasing damage. However, both of these are 
macroscopic quantities. The microscopic level of the process is of greater in- 
terest if we want to understand how the fragmentation occurs. In our model, 
cracks and fragments are formed by the breaking of bonds between neighboring 
particles. This is illustrated in the snapshots of Fig. where the bonds that 
still exist at the end of the fragmentation process (t = 100) are shown in their 
initial locations before the impact. Fragments can be discerned as clusters of 
the remaining bonds (except the monomers, of course). 

At low initial kinetic energy (V — 1, i] — 0.17), the disc suffers mainly 
from internal damage. Some internal bonds are broken and a few fragments are 
formed, mainly monomers in the impact zone (bottom of the figure). Oblique 
lines of broken bonds starting from the impact zone can be seen. As the energy 
deposited in the system increases (V = 2, r\ = 0.70), oblique cracks propagate 
through the system from the impact zone to the edge of the disc and fragmen- 
tation in the proper sense of the word occurs. The largest fragment is localized 
in the top part of the disc. Most of the energy is dissipated in the impact zone, 
producing quantities of monomers. At higher energies (V > 3 , r\ > 1.57), 
cracks perpendicular to the former oblique ones form and create smaller frag- 
ments by merging. As energy increases further, fragments become smaller but 
are of various sizes, in qualitative agreement with a power law fragment size dis- 
tribution. It is interesting to note, that the shape of the fragments is irregular 
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Figure 3: Snapshots at the end of the fragmentation process (t — 100 ) for various 
impact velocities and the same initial rotation angle. Each dot represents a 
bond between two particles localized at its initial position before the impact. The 
horizontal wall is depicted at the bottom of the top left figure. 



and rough without memory of the underlying triangular lattice structure of the 
disc. To conclude, the crack patterns shown in Fig. [3] are consistent with the 
ones obtained with more complex models [TBI [2Tj . 

To quantify the damage produced during the impact fragmentation, Thorn- 
ton and his collaborators [TH] defined the damage ratio D(t) as the average ratio 
of the number of broken bonds at a given time t, to the total number of initial 
bonds in the disc. This definition of damage is quite natural at a microscopic 
scale and can be easily applied to many kinds of fragmenting systems. In the 
inset of Fig. 2] we plot the damage ratio for various impact velocities. At the 
moment of impact (t = 0), D{t) first increases dramatically, then at a slower 
rate, and eventually tends to its final value Df, showing that the fragmentation 
process reaches a steady state. In fact, at later times (t > 50), only few particles 
evaporate from exited fragments, increasing very slightly the damage ratio. 

In order to better understand the evolution of damage, we calculated the 

damage rate, ^ . As shown in Fig 01 the damage rate increases drastically, 
dt 

reaches a maximum, and then decreases slowly to zero. 

This evolution is further illustrated by the snapshots in Fig [5] of the frag- 
menting disc for the impact velocity V — 3 (n — 1.57). Just after the impact, 
at t\ = 0.15, two oblique cracks start to propagate from the impact zone. At 
t% = 0.35, damage increases in the impact zone and the oblique cracks become 
wider and continue their advance through the disc. As a consequence the dam- 
age rate increases swiftly. When the two cracks have reached the edge of the 
disc, the damage rate is maximal. Then, cracks perpendicular to the oblique 
ones are formed {t$ = 0.55). Damage spreads into the disc and damage rate 
decreases (t > £3). The local minimum of the damage rate (£2 < t < £3) is due 
to a sudden increase of the kinetic energy per particle (not shown) when the 
disc rebounds, after being compressed (this is observed also for higher impact 
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Figure 4: Damage rate ) as a function of time for various impact velocities. 
The times U, for i = 1,2,3,4,5, along the curve for V = 3 correspond to the 
snapshots of Fig. In the inset, D(t) as a function of time for the same impact 
velocities. 



velocities, see the curve lor V = 5). Finally, bond breakings create fragments 
whose surfaces evolve by evaporating particles (t > t$ = 1.55). As can be seen 
by comparing Fig. [5] at f 6 = 10 and Fig. [3] at t — 100, fragments already have 
their overall shape at t = 10. Note that the lattice structure influences the first 
oblique crack propagation. However, these linear cracks are not sufficient to 
dissipate the impact energy and secondary winding cracks form within the disc 
without revealing the underlying lattice structure. 

As illustrated by the curves in Fig. [H as the impact velocity increases, the 
damage rate increases faster and reaches a maximal value when the oblique 
cracks have propagated through the system. The position of the peak depends 
only slightly on the energy deposited in the system, while it is mostly determined 
by the direction of propagation of the cracks. Indeed, as the speed of sound in 
the disc is estimated to be around 100 [53], a crack of length 40cr (the diameter 
of the disc) reaches the edge at time t ~ 0.4. This is in agreement with the 
position of the peaks in Fig. 0] Note that the regime of evaporation is reached 
at the same time, around t ~ 1, regardless of the energy deposited in the system. 

Having studied the evolution of damage, we discuss how the final damage 
ratio Df at the end of the process behaves as a function of the initial kinetic 
energy. As illustrated by the curve in Fig. the final damage ratio Df is 
very well fitted by a logarithmic function for 0.4 < n < 6, or equivalently for 
0.2 <D f < 0.9: 

D f = a Df bx^-. (4) 

Here au } — 0.256 ± 0.005, i.e., almost the same coefficient as for r and for 
mo/N. For higher energies (r) > 6), Df grows more slowly than the logarithmic 
function, reflecting a lower efficiency of the fragmentation process. Therefore, 
even at very high initial kinetic energies, some small fragments still remain. For 
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Figure 5: Snapshots at various times for V — 3 (rj — 1.57). Each dot represents 
a bond between two particles placed at its initial position before the impact. 

instance, at V — 8 (n = 11.15), S max i — 13. 

The same logarithmic behaviour of Df was reported by Thornton et al. [18] 
and Bchera et al. [21j . For the former model, they found ao f — 0.1, depending 
on the surface energy of the disc. For the latter one, otr> s — 0.15, as can be 
estimated from Fig. 8 of reference [2"T] . 

According to equations J2J), ©, and the power law exponent r, the total 
number of fragments mo/N, and the damage ratio Df all exhibit a logarithmic 
dependence on r\ with the same coefficient a, in a certain range of energy. 
Consequently, we expect a linear relation between r and Df, and between mo/N 
and Df. This is plotted in the inset of Fig. [BJ A linear fit gives 

r = D/ + 1.35, (5) 

for 0.2 < Df < 0.7 and 

m /N = 0.961?/ -0.16, (6) 

for 0.2 < Df < 0.9, with a slope slightly less than 1. Note that using the data 
presented by Behera et al. |21j , we found a linear relation between mo and D / 
as well. 



5 Fragmentation energy 

A large part of the total energy consumed in industry is used to achieve size 
reduction of materials. To estimate the efficiency of the fragmentation process 
studied in this article, we calculated the fragmentation energy. Noting that no 
frictional interaction is taken into account in the present model, it follows from 
the conservation of energy that 

1 2 

e po ta + i^mV = e pot + e kin , (7) 
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Figure 6: Final damage ratio Df at t = 100 as a function of rj. The full line is 
a logarithmic fit (see eq. RW. The power law exponent t and the total number 
of fragments mo/N are plotted in the inset as functions of damage ratio Df (the 
full lines are linear fits). 



where e po t and ekin are the average potential and kinetic energies per particle at 
the end of the fragmentation process, and e pot o = — 2.87e is the initial potential 
energy per particle. In other words, the initial kinetic energy is used to break 
bonds and create fragments (potential energy) and to move the fragments apart 
(kinetic energy). Consequently, we define the fragmentation energy ef rag as the 
potential energy needed to form the fragments: ef rag = e pot — e po to- 

It is easy to see that the fragmentation energy is related to the final damage 
ratio Df. Using a linear fit, we found for Df > 0.1: 

e frag = 2.75D f + 0.12. (8) 

Due to surface effects, the fragmentation energy is lower than —e po toDf when 
the disc is broken into fragments. When the disc is completely shattered and 
Df = 1, we have ef rag = —e po to, as expected. 

The efficiency of the fragmentation process can be estimated by calculating 
the fraction of energy used to create the fragments, that is, e f rag / ekino where 
efcino = \mV 2 is the initial kinetic energy per particle. This quantity is plotted 
in Fig. [7] as a function of r) in the fragmentation regime, i.e., for r\ > rj t - The 
percentage of the initial kinetic energy that is used to form the fragments de- 
creases from 70% to 10% as r\ increases from r\ t to 10. At high impact velocities, 
most of the energy is then spent on moving the fragments apart. If friction 
would be taken into account, the efficiency of the fragmentation process would 
be even lower. 
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Fi gure 7: Fragmentation energy divided by the initial kinetic energy, ef rag / e^mOi 
as a function of r), for rj > r/ t . 

6 Conclusion 

The model we propose is generic in the sense that it contains the essential phys- 
ical features for investigating fragmentation phenomena. Its simplicity rests on 
a basic cohesive interaction between circular particles in two dimensions. Fur- 
thermore, molecular dynamics calculations allow us to determine the physical 
quantities involved and moreover, to study their evolution during the fragmen- 
tation process. 

In this article, we have investigated the damage inflicted on a disc of inter- 
acting particles as the disc fragments upon impact with a wall — the damage 
being defined as the percentage of broken bonds between particles. Its propaga- 
tion into the disc evidently depends on the energy deposited in the system, but 
the damage rate was found to be independent of the impact energy shortly after 
the collision. Furthermore, the percentage of energy actually used to fragment 
the system decreases with the energy deposited in the system, reflecting the low 
efficiency of the fragmentation process at high impact energy. 

In particular, we have shown that the power law exponent t of the fragment 
size distribution, the total number of fragments, and the damage share the 
same logarithmic behaviour as functions of the energy deposited in the system. 
Damage cannot be directly measured by experiments, but the logarithmic be- 
haviour of r is observed in the fragmentation of rocks [2] and liquid droplets [10] . 
Surprisingly, these two very different systems exhibit the same fragmentation 
behaviour. This suggests that at high impact energies, the crack propagation 
and fragmentation phenomena may not be sensitive to the particular structure 
of the system. As an interesting consequence of this same logarithmic behaviour, 
the power law exponent and the total number of fragments are proportional to 
the damage. These simple relations, that we expect to be generic, may guide 
us to a better understanding of the underlying mechanisms of fragmentation 
phenomena. 
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